home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Tech Arsenal 1
/
Tech Arsenal (Arsenal Computer).ISO
/
tek-20
/
ffteme12.zip
/
FFTEME12.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1991-04-02
|
11KB
|
498 lines
{$R+} {Range checking}
{$B+} {Boolean complete evaluation}
{$S+} {Stack checking}
{$I+} {I/O checking}
{$N+} {numeric coprocessor}
program FFTEME;
{ by Mike Cook (AF9Y) Mar 31, 1991}
Uses
DOS,
Crt,
Graph;
var
GraphDriver: integer;
GraphMode: integer;
ErrorCode: integer;
FileVar: File of Byte;
FastFile: File;
OutFile: Text;
FileName: String[14];
SampleByte: Byte;
ImagF: Array [0..4096] of Single;
RealF: Array [0..4096] of Single;
Sample: Array [0..28000] of Byte;
Mag: Single;
Period: Single;
TotalSamples: Integer;
TotalFileSamples: Integer;
Ch: Char;
N: Integer;
NLimit: Integer;
X1: Integer; X2: Integer;
Y1: Integer; Y2: Integer;
I,X,L,Y,K,M,MX: Integer;
LoopEnd: Integer;
X3,Max: Single;
Z,Zmax: Integer;
ADwait: Integer;
GraphLine: Integer;
PreviousGraphLine: Integer;
YScale: Single;
YStart: Integer;
XOffSet: Integer;
YOffSet: Integer;
NormFFTAns: Char;
NormFFT: Integer;
Start: Integer;
StringZ: String[40];
DefaultMinFreq, MinFreq: Integer;
DefaultMaxFreq, MaxFreq: Integer;
MinZ: Integer;
MaxZ: Integer;
DataSkips: Integer;
Skips: Integer;
DisplayRow: Integer;
DisplayCol: Integer;
HourA,MinuteA,SecondA,Sec100A: Word;
HourB,MinuteB,SecondB,Sec100B: Word;
StartTime, EndTime: LongInt;
Break: Boolean;
Procedure GraphSetup;
Begin
GraphDriver:=Detect;
InitGraph(GraphDriver, GraphMode, '');
end;
Procedure InputSampleTotal;
Var
OffsetLimit: Integer;
Begin
Write ('Number of Samples (4096,2048,1024)? '); Readln(TotalSamples);
L:=Round(Ln(TotalSamples)/Ln(2));
Period:=((TotalSamples/4096)*0.2485);{Recorded Sample Rate Cal aginist wwv}
ADwait:=Round(6+(28500*Period/TotalSamples));
Write ('DataSkips (1,2,3,4)? '); Readln(DataSkips);
Period:=Period*DataSkips;
OffsetLimit:=Trunc(27776-(DataSkips*TotalSamples));
Write ('Offset Samples (0..',OffsetLimit,')? '); Readln(Start);
Write ('Normal or Squared FFT (N/S)? '); Readln(NormFFTAns);
If ((NormFFTAns='S') or (NormFFTAns='s')) then NormFFT:=0 else NormFFT:=1;
If NormFFT=0 then Period:=Period*2;
DefaultMinFreq:=Round(1/Period);
DefaultMaxFreq:=Round(TotalSamples/2/Period);
WriteLn('FFT Range = ',DefaultMinFreq,' Hz',' to ',DefaultMaxFreq,' Hz');
Write ('Min FFT Freq? '); Readln(MinFreq);
MinZ:=Round(MinFreq*Period);
Write ('Max FFT Freq? '); Readln(MaxFreq);
MaxZ:=Round(MaxFreq*Period);
End;
Procedure RequestReadFile;
begin
Write('Enter name of file to read: ');
Readln(FileName);
{FileName:='EME1';}
Assign(FileVar,FileName);
Reset(FileVar);
end;
Procedure FastReadFile;
Var
NumRead: Word;
Begin
Assign(FastFile,FileName);
Reset(FastFile,1);
Repeat
BlockRead(FastFile,Sample,SizeOf(Sample),NumRead);
Until (NumRead = 0);
Close(FileVar);
NLimit:=TotalSamples-1;
end;
Procedure ReadFile;
Var
ReadCnt: Integer;
TotalReads: Single;
StoreCnt: Integer;
Begin
WriteLn;
DisplayRow:=WhereY;
DisplayCol:=WhereX;
StoreCnt:=0; ReadCnt:=0; TotalReads:=(TotalSamples*DataSkips+Start)/100;
While ((Not EOF(FileVar)) AND (StoreCnt<TotalSamples)) do
Begin
If ReadCnt<Start then
Begin
Read(FileVar,SampleByte);
Inc(ReadCnt);
ClrEol; GotoXY(DisplayCol,DisplayRow);
Write('Data Read ',(Round(ReadCnt/TotalReads)),'% Complete');
end
else begin
For Skips:=1 to DataSkips do
Begin
Read(FileVar,SampleByte);
Inc(ReadCnt);
If EOF(FileVar) then Skips:=DataSkips;
end;
Sample[(StoreCnt)]:=(SampleByte);
StoreCnt:=StoreCnt+1;
ClrEol; GotoXY(DisplayCol,DisplayRow);
Write('Data Read ',(Round(ReadCnt/TotalReads)),'% Complete');
end;
end;
Close(FileVar);
If StoreCnt<TotalSamples then
Begin
Writeln('Not Enough Data');
Halt;
end;
End;
Procedure DisableClock;
Begin
Port[$21]:=(Port[$21] or $01);
End;
Procedure EnableClock;
Begin
Port[$21]:=(Port[$21] and $FE);
End;
Procedure GenSineWave;
Var
Freq: Single;
TimeStep: Single;
TimeN: Single;
Period: Single;
Begin
Freq:=555; TimeStep:=Period/Totalsamples;
TimeN:=0;
NLimit:=TotalSamples-1;
For N:=0 to Nlimit do
Begin
Sample[N]:=Round(128*(Sin(2*PI*Freq*TimeN)));
TimeN:=TimeN+TimeStep;
End;
End;
Procedure CopyFastSamples;
Var
IndexN: Integer;
Begin
NLimit:=TotalSamples-1;
IndexN:=Start;
For N:=0 to NLimit do
Begin
If NormFFT=1 then RealF[N]:=((Sample[IndexN]-128)/128)/TotalSamples
{ Normal FFT }
else begin
RealF[N]:=(((Sample[IndexN]-128)*(Sample[IndexN]-128))/128-128)/TotalSamples;
{ Square FFT }
end;
ImagF[N]:=0;
IndexN:=IndexN+DataSkips;
end;
End;
Procedure CopySamples;
Begin
NLimit:=TotalSamples-1;
For N:=0 to NLimit do
Begin
If NormFFT=1 then RealF[N]:=(Sample[N]/128-128)/TotalSamples { Normal FFT }
else
RealF[N]:=((Sample[N])*Sample[N]/128-128)/TotalSamples; { Square FFT }
ImagF[N]:=0;
end;
End;
Procedure Scramble;
Var
W: Integer;
N1: Integer;
Begin
Y:=0; N1:=TotalSamples;
For W:=1 to L do
Begin
N1:=N1 shr 1;
If X>=N1 then
Begin
If W=1 then Inc(Y)
else
Y:=Y+(2 shl (W-2));
X:=X-N1;
end;
end;
End;
Procedure CalFFT;
Var
I1,I2,I3,I4,I5: Integer;
V,A1,A2,B1,B2,Z1,Z2: Single;
Begin
I1:=TotalSamples shr 1; I2:=1; V:=2*Pi/TotalSamples;
WriteLn;
DisplayRow:=WhereY;
DisplayCol:=WhereX;
For I:=1 to L do
Begin
ClrEol; GotoXY(DisplayCol,DisplayRow);
Write('FFT ',(Round((I/L)*100)),'% Complete');
I3:=0; I4:=I1;
For K:=1 to I2 do
Begin
X:=Trunc(I3/I1); Scramble;
I5:=Y;
Z1:=Cos(V*I5); Z2:=-Sin(V*I5);
LoopEnd:=I4-1;
For M:=I3 to LoopEnd do
Begin
A1:=RealF[M]; A2:=ImagF[M]; Mx:=M+I1;
B1:=Z1*RealF[Mx]-Z2*ImagF[Mx];
B2:=Z2*RealF[Mx]+Z1*ImagF[Mx];
RealF[M]:=(A1+B1); ImagF[M]:=(A2+B2);
RealF[Mx]:=(A1-B1); ImagF[Mx]:=(A2-B2);
end;
I3:=I3+(I1 shl 1); I4:=I4+(I1 shl 1);
end;
I1:=I1 shr 1; I2:=I2 shl 1;
end;
end;
Procedure CalMagnitude;
Begin
Scramble;
Mag:=SQRT(RealF[Y]*RealF[Y]+ImagF[Y]*ImagF[Y]);
End;
Procedure CalMaxSpectrum;
Begin
Max:=0;
LoopEnd:=TotalSamples shr 1;
For Z:=MinZ to MaxZ do
Begin
X:=Z;
CalMagnitude;
If Mag>Max then
Begin
Max:=Mag;
Zmax:=Z;
End;
End;
End;
Procedure PlotTimeSamples;
Begin
YScale:=0.13; YStart:=10;
GraphLine:=0; PreviousGraphLine:=0;
LoopEnd:=TotalSamples-1;
For N:=0 to LoopEnd do
Begin
Case N of
0000..0639 : GraphLine:=1;
0640..1279 : GraphLine:=2;
1280..1919 : GraphLine:=3;
1920..2559 : GraphLine:=4;
2560..3199 : GraphLine:=5;
3200..3839 : GraphLine:=6;
3840..4479 : GraphLine:=7;
4480..5120 : GraphLine:=8;
end;
If Graphline<>PreviousGraphLine then
Begin
YOffSet:=YStart+GraphLine*35;
Y1:=Round(YOffSet-YScale*(Sample[N]));
XOffSet:=Round((GraphLine-1)*640);
PreviousGraphLine:=GraphLine;
end;
Y2:=Round(YOffSet-YScale*(Sample[N]+128));
Line(N-XOffSet,Y1,N-XOffSet,Y2);
Y1:=Y2
End;
end;
Procedure PlotFastTimeSamples;
Var
IndexN: Integer;
Begin
YScale:=0.13; YStart:=10;
GraphLine:=0; PreviousGraphLine:=0;
LoopEnd:=TotalSamples-1;
For N:=0 to LoopEnd do
Begin
Case N of
0000..0639 : GraphLine:=1;
0640..1279 : GraphLine:=2;
1280..1919 : GraphLine:=3;
1920..2559 : GraphLine:=4;
2560..3199 : GraphLine:=5;
3200..3839 : GraphLine:=6;
3840..4479 : GraphLine:=7;
4480..5120 : GraphLine:=8;
end;
IndexN:=N*DataSkips+Start;
If Graphline<>PreviousGraphLine then
Begin
YOffSet:=YStart+GraphLine*35;
Y1:=Round(YOffSet-YScale*(Sample[IndexN]));
XOffSet:=Round((GraphLine-1)*640);
PreviousGraphLine:=GraphLine;
end;
Y2:=Round(YOffSet-YScale*(Sample[IndexN]+128));
Line(N-XOffSet,Y1,N-XOffSet,Y2);
Y1:=Y2
End;
end;
Procedure PlotSpectrum;
Begin
YScale:=0.40; YStart:=260;
GraphLine:=0; PreviousGraphLine:=0;
LoopEnd:=TotalSamples Shr 1;
For Z:=MinZ to MaxZ do
Begin
Case Z of
0000..0639 : GraphLine:=1;
0640..1279 : GraphLine:=2;
1280..1919 : GraphLine:=3;
1920..2559 : GraphLine:=4;
2560..3200 : GraphLine:=5;
3200..3839 : GraphLine:=6;
3840..4480 : GraphLine:=7;
4480..5120 : GraphLine:=8;
end;
If Graphline<>PreviousGraphLine then
Begin
YOffSet:=Round(YStart+GraphLine*50.0);
XOffSet:=Round((GraphLine-1)*640.0);
PreviousGraphLine:=GraphLine;
end;
X:=Z;
CalMagnitude;
Y1:=Round(YOffSet-(YScale*Mag/Max)*199.0);
Line(Z-XOffSet,YOffSet,Z-XOffSet,Y1);
End;
end;
Procedure PrintSpectrum;
Begin
writeln;
LoopEnd:=TotalSamples shr 1;
For Z:=0 to LoopEnd do
Begin
X:=Z;
CalMagnitude;
WriteLn(Z,' ',RealF[Y],' ',ImagF[Y],' ',Mag);
End;
end;
{============================= MAIN ====================================}
Begin
Break:=False;
Repeat
Begin
TextMode(CO80);
ClrScr;
InputSampleTotal;
{GenSineWave;}
RequestReadFile;
{ReadFile;}
FastReadFile;
{CopySamples;}
CopyFastSamples;
GetTime(HourA,MinuteA,SecondA,Sec100A);
CalFFT;
GetTime(HourB,MinuteB,SecondB,Sec100B);
StartTime:=HourA*360000+MinuteA*6000+SecondA*100+Sec100A;
EndTime:=HourB*360000+MinuteB*6000+SecondB*100+Sec100B;
Writeln;
Writeln('FFT Process Time = ',((EndTime-StartTime)/100):2:2);
CalMaxSpectrum;
GraphSetup;
PlotFastTimeSamples;
{PlotTimeSamples;}
{PrintSpectrum;}
PlotSpectrum;
Zmax:=Round(Zmax/Period);
Str(Zmax,StringZ); OutTextXY(600,460,StringZ);
Repeat Until Keypressed;
end;
Until Break;
TextMode(CO80);
End.